This document includes the statistical analyses reported in the paper Dablander, F.â‘, Sachisthal, M.â‘, & Aron, A. (under review) Out of the Labs and into the Streets? Effects of Climate Protests by Environmental Scientists. Royal Society Open Science.
We prepare the data and then test the above hypotheses.
library(brms)
library(dplyr)
library(knitr)
library(tidyr)
library(readr)
library(scales)
library(ggplot2)
library(corrplot)
library(gridExtra)
library(BayesFactor)
library(marginaleffects)
source('helpers.R')
df_text <- read_csv('data/dat_text.csv')
df <- read_csv('data/dat_num_all.csv')
demographics <- read.csv('data/demographics_all.csv')
df <- df %>%
rename(Submission.id = SESSION_ID) %>%
left_join(
demographics %>%
rename(political_affiliation = U.s..political.affiliation) %>%
select(Submission.id, Age, Sex, political_affiliation),
by = 'Submission.id'
) %>%
mutate(
Protest = case_when(
grepl('Legal', condition) ~ 'Legal March',
grepl('CD', condition) ~ 'Civil Disobedience'
),
Protest = factor(Protest, levels = c('Civil Disobedience', 'Legal March')),
Engagement = case_when(
grepl('None', condition) ~ 'Control',
grepl('Endorse', condition) ~ 'Support',
grepl('Join', condition) ~ 'Engaged'
),
# Engagement = factor(Engagement, levels = c('Engaged', 'Support', 'Control')),
Engagement = factor(Engagement, levels = c('Control', 'Support', 'Engaged')),
PolicySupport = as.numeric(PolicySupport),
ActivistSupport = as.numeric(ActivistSupport),
ScienceCredibility = as.numeric(ScienceCredibility),
SourceCredibility = as.numeric(SourceCredibility),
Radical = as.numeric(Radical),
Donation = as.numeric(Donation_1),
PoliticalAffiliation = as.numeric(PoliticalAffiliation)
) %>%
filter(
# !(political_affiliation %in% c('CONSENT_REVOKED', 'DATA_EXPIRED')),
Progress == 100 # Remove 10 people who did not finish the survey
) %>%
mutate(
political_affiliation = ifelse(
PoliticalAffiliation %in% c(1, 2), 'Democrat',
ifelse(
PoliticalAffiliation %in% c(4, 5), 'Republican', 'Independent'
)
),
Politics = factor(political_affiliation, levels = c('Democrat', 'Republican', 'Independent'))
)
contrasts(df$Protest) <- contr.sum(levels(df$Protest))
contrasts(df$Engagement) <- contr.sum(levels(df$Engagement))
contrasts(df$Politics) <- contr.sum(levels(df$Politics))
nrow(df)
## [1] 3149
We remove those participants that incorrectly answer the attention check question.
df$AC_Protest %>% table
## .
## 1 2 3
## 1755 1350 44
df_text$AC_Protest %>% table
## .
## Blockade with arrests Peaceful protest without arrests
## 1350 1754
## There was no protest
## 44
df <- df %>%
filter(
(Protest == 'Civil Disobedience' & AC_Protest == '2') |
(Protest == 'Legal March' & AC_Protest == '1')
)
We are left with 2,875 responses.
nrow(df)
## [1] 2875
We remove the users that responded incoherently or inappropriately to the open text question about what the article was about. This leaves us with 2,856 responses.
df_open <- read.csv('data/openAnswers_coded.csv') %>% filter(Usable == 1)
# Additional check from second coder
df_check <- df_open %>% filter(Check == 1)
remove_responseids <- df_check$ResponseId[-c(4, 6)]
df <- df %>%
left_join(df_open, by = 'ResponseId') %>%
filter(Usable == 1, !(ResponseId %in% remove_responseids))
nrow(df)
## [1] 2856
# Prepare sub-set of data for later analyses
df_march <- df %>%
filter(Protest == 'Legal March') %>%
mutate(Y = Radical)
df_civil <- df %>%
filter(Protest == 'Civil Disobedience') %>%
mutate(Y = Radical)
df_source <- df %>%
filter(!is.na(SourceCredibility)) %>%
mutate(
Y = SourceCredibility,
Engagement = factor(Engagement, levels = c('Support', 'Engaged'))
)
contrasts(df_source$Engagement) <- contr.sum(levels(df_source$Engagement))
We fit all the models for the model-averaged exploratory analysis.
# Not included because it is ~50 MB large
fit_all <- readRDS('models/models_scaled.RDS')
# Not included because it is ~300 MB large
if (file.exists('models/fit_all_averaging_scaled.RDS')) {
# For model-averaging later
fit_all_averaging <- readRDS('models/fit_all_averaging_scaled.RDS')
fit_policy <- fit_all_averaging$fit_policy
fit_activist <- fit_all_averaging$fit_activist
fit_radical <- fit_all_averaging$fit_radical
fit_source <- fit_all_averaging$fit_source
fit_science <- fit_all_averaging$fit_science
} else {
set.seed(1)
fit_policy <- get_logml_all_parallel(
df %>% mutate(Y = PolicySupport), fit_all,
iter = 4000, cores = 1, chains = 2
)
set.seed(1)
fit_activist <- get_logml_all_parallel(
df %>% mutate(Y = ActivistSupport), fit_all,
iter = 4000, cores = 1, chains = 2
)
set.seed(1)
fit_radical <- get_logml_all_parallel(
df %>% mutate(Y = Radical), fit_all,
iter = 4000, cores = 1, chains = 2
)
set.seed(1)
fit_source <- get_logml_all_parallel(
df_source, fit_all,
iter = 4000, cores = 1, chains = 2, recompile = TRUE
)
set.seed(1)
fit_science <- get_logml_all_parallel(
df %>% mutate(Y = ScienceCredibility), fit_all,
iter = 4000, cores = 1, chains = 2
)
fit_all_averaging <- list(
'fit_policy' = fit_policy,
'fit_activist' = fit_activist,
'fit_radical' = fit_radical,
'fit_source' = fit_source,
'fit_science' = fit_science
)
saveRDS(fit_all_averaging, 'models/fit_all_averaging_scaled.RDS')
}
fit <- fit_policy$fit_updated[[5]]
logmls_policy <- fit_policy$logmls
df_effects_policy <- get_all_predictions(
fit, 'Policy support', variables = c('Engagement', 'Protest')
)
palette <- c("#E67C73", "#F2A489", "#F6E8A7", "#A3C4D4", "#5D92C7")
p1 <- create_plot(
df, 'PolicySupport', 'Policy support',
df_effects = df_effects_policy
)
p1
Moderate evidence against effect of protest.
h1_bfs <- get_bf_matrix(logmls_policy[c(1, 3)])
h1_bfs
## H0 H1 H1r
## H0 vs 1.0000000 8.05561 NA
## H1 vs 0.1241371 1.00000 NA
## H1r vs NA NA NA
Strong evidence of no effect of engagement on policy support.
h4_bfs <- get_bf_matrix(
c(logmls_policy[c(1, 2)],
log(get_directed_bf(fit_policy$fit_updated[[2]])) + logmls_policy[2]
)
)
h4_bfs
## H0 H1 H1r
## H0 vs 1.00000000 38.89893 43.953589
## H1 vs 0.02570765 1.00000 1.129944
## H1r vs 0.02275127 0.88500 1.000000
fit <- fit_activist$fit_updated[[5]]
logmls_activist <- fit_activist$logmls
df_effects_activist <- get_all_predictions(
fit, 'Activist support', variables = c('Engagement', 'Protest')
)
p2 <- create_plot(
df, 'ActivistSupport',
'Activist support',
df_effects = df_effects_activist
)
p2
Overwhelming evidence in favour of a difference.
h2_bfs <- get_bf_matrix(
c(logmls_activist[c(1, 3)],
logmls_activist[3] + log(get_directed_bf(
fit_activist$fit_updated[[3]], groups = 2,
type = 'decreasing', pars = 'Protest'
))
)
)
h2_bfs
## H0 H1 H1r
## H0 vs 1.000000e+00 1.745622e-15 8.728109e-16
## H1 vs 5.728617e+14 1.000000e+00 5.000000e-01
## H1r vs 1.145723e+15 2.000000e+00 1.000000e+00
Strong evidence in favour of no effect, although it’s a bit lower compared to the directed hypothesis.
h5_bfs <- get_bf_matrix(
c(logmls_activist[c(1, 2)],
logmls_activist[2] + log(get_directed_bf(
fit_activist$fit_updated[[2]]
))
)
)
h5_bfs
## H0 H1 H1r
## H0 vs 1.00000000 14.43159 4.297033
## H1 vs 0.06929246 1.00000 0.297752
## H1r vs 0.23271871 3.35850 1.000000
fit <- fit_radical$fit_updated[[5]]
logmls_radical <- fit_radical$logmls
df_effects_radical <- get_all_predictions(
fit, 'Perceived radicalness', variables = c('Engagement', 'Protest')
)
p3 <- create_plot(
df, 'Radical',
'Perceived radicalness',
df_effects = df_effects_radical,
lo = 'Not at all radical', hi = 'Extremely radical'
)
p3
Overwhelming effect of higher perceived radicalness for civil disobedience.
h3_bfs <- get_bf_matrix(
c(logmls_radical[c(1, 3)],
logmls_radical[3] + log(get_directed_bf(
fit_radical$fit_updated[[3]], groups = 2,
type = 'increasing', pars = 'Protest'
))
)
)
h3_bfs
## H0 H1 H1r
## H0 vs 1.000000e+00 6.036804e-79 3.018402e-79
## H1 vs 1.656506e+78 1.000000e+00 5.000000e-01
## H1r vs 3.313011e+78 2.000000e+00 1.000000e+00
Strong evidence evidence in favour of no effect of engagement.
h6_test <- get_logml_all(
df_march,
list('null' = fit_all[[1]], 'full' = fit_all[[2]])
)
## [1] 0.15
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
# All marginal likelihoods
logmls_all <- c(
h6_test$logmls,
h6_test$logmls[2] + log(get_directed_bf(h6_test$fit_updated$full, type = 'decreasing'))
)
h6_bfs <- get_bf_matrix(logmls_all)
h6_bfs
## H0 H1 H1r
## H0 vs 1.00000000 4.301729 37.73446
## H1 vs 0.23246468 1.000000 8.77193
## H1r vs 0.02650097 0.114000 1.00000
Weak evidence evidence in favour of no effect of engagement.
h7_test <- get_logml_all(
df_civil,
list('null' = fit_all[[1]], 'full' = fit_all[[2]])
)
## [1] 0.15
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6
# All marginal likelihoods
logmls_all <- c(
h7_test$logmls,
h7_test$logmls[2] + log(get_directed_bf(h7_test$fit_updated$full, type = 'decreasing'))
)
h7_bfs <- get_bf_matrix(logmls_all)
h7_bfs
## H0 H1 H1r
## H0 vs 1.0000000 0.6879398 1.777622
## H1 vs 1.4536157 1.0000000 2.583979
## H1r vs 0.5625493 0.3870000 1.000000
We need to recompile the model because otherwise
marginaleffects throws an error regarding the different
factor levels for engagement being present.
# Is small enough to be included, but also quick enough to be estimated once
if (!file.exists('models/fit_source_full_final.RDS')) {
fit <- brm(
formula = 'Y ~ Engagement + Protest',
data = df_source, family = brms::cumulative('probit'),
prior = c(
prior(student_t(1, 0, 0.15), class = b),
prior(student_t(1, 0, 2.5), class = Intercept)
), cores = 4, chains = 4, iter = 1000, save_pars = save_pars(all = TRUE)
)
saveRDS(fit, 'models/fit_source_full_final.RDS')
} else {
fit <- readRDS('models/fit_source_full_final.RDS')
}
df_effects_source <- get_all_predictions(
fit, 'Source credibility', variables = c('Engagement', 'Protest')
) %>%
filter(Engagement != 'Control')
p4 <- create_plot(
df_source, 'SourceCredibility',
'Source credibility',
df_effects = df_effects_source,
lo = 'Not at all', hi = 'A great deal'
)
p4
Weak evidence in favour of no effect.
logmls_source <- fit_source$logmls
h8_bfs <- get_bf_matrix(logmls_source[c(1, 2)])
h8_bfs
## H0 H1 H1r
## H0 vs 1.0000000 1.854321 NA
## H1 vs 0.5392809 1.000000 NA
## H1r vs NA NA NA
fit <- fit_science$fit_updated[[5]]
df_effects_science <- get_all_predictions(
fit, 'Science credibility', variables = c('Engagement', 'Protest')
)
p5 <- create_plot(
df, 'ScienceCredibility',
'Science credibility',
df_effects = df_effects_science,
lo = 'Not at all', hi = 'A great deal'
)
p5
Strong evidence for no effect of engagement.
logmls_science <- fit_science$logmls
h9_bfs <- get_bf_matrix(logmls_science[c(1, 2)])
h9_bfs
## H0 H1 H1r
## H0 vs 1.00000000 38.4238 NA
## H1 vs 0.02602553 1.0000 NA
## H1r vs NA NA NA
ggsave(
'figures/Figure1.pdf',
grid.arrange(
p1, p2, p3, p5, p4,
layout_matrix = rbind(
c(1, 1, 2, 2),
c(3, 3, 4, 4),
c(NA, 5, 5, NA)
)
),
width = 11, height = 12
)
hyp_names <- c(
'H1: No effect of protest form on policy support',
'H2: Higher activist support for legal march',
'H3: Higher perceived radicalness for civil disobedience',
'H4: Higher policy support for scientist joining, then supporting, then for no engagement',
'H5: Higher activist support for scientist joining, then supporting, then for no engagement',
'H6: Legal march perceived as least radical when scientists join, then for supporting, then for no engagement',
'H7: Civil disobedience perceived as least radical when scientists join, then for supporting, then for no engagement',
'H8: No effect of type of engagement on source credibility',
'H9: No effect of type of engagement on science credibility'
)
hyp_bfs <- c(
h1_bfs[1, 2],
h2_bfs[3, 1],
h3_bfs[3, 1],
h4_bfs[3, 1],
h5_bfs[3, 1],
h6_bfs[3, 1],
h7_bfs[3, 1],
h8_bfs[1, 2],
h9_bfs[1, 2]
)
kable(
cbind(hyp_names, round(hyp_bfs, 3)),
col.names = c('Hypothesis', 'Bayes factor in favour of the hypothesis')
)
| Hypothesis | Bayes factor in favour of the hypothesis |
|---|---|
| H1: No effect of protest form on policy support | 8.056 |
| H2: Higher activist support for legal march | 1145723482010399 |
| H3: Higher perceived radicalness for civil disobedience | 3.31301116447911e+78 |
| H4: Higher policy support for scientist joining, then supporting, then for no engagement | 0.023 |
| H5: Higher activist support for scientist joining, then supporting, then for no engagement | 0.233 |
| H6: Legal march perceived as least radical when scientists join, then for supporting, then for no engagement | 0.027 |
| H7: Civil disobedience perceived as least radical when scientists join, then for supporting, then for no engagement | 0.563 |
| H8: No effect of type of engagement on source credibility | 1.854 |
| H9: No effect of type of engagement on science credibility | 38.424 |
Here we show the results across political affiliation and report inclusion Bayes factors.
df_effects_policy <- get_all_predictions(
fit_policy$fit_updated[[19]], 'Policy support',
variables = c('Engagement', 'Protest', 'Politics')
)
p7 <- create_plot(
df, 'PolicySupport',
'Policy support',
'Expand offshore drilling',
political = TRUE,
df_effects = df_effects_policy
)
p7
ggsave('figures/FigureS4.pdf', p7, width = 9, height = 7)
inclusion_bfs_policy_support <- get_inclusion_bayes_factors(fit_policy)
kable(
log(inclusion_bfs_policy_support),
col.names = c('Effect', 'Log Inclusion Bayes factor')
)
| Effect | Log Inclusion Bayes factor |
|---|---|
| Engagement | -3.3143482 |
| Protest | -0.9983041 |
| Politics | 365.1042998 |
| Engagement:Protest | -1.1661951 |
| Engagement:Politics | -2.8801898 |
| Protest:Politics | -2.2245190 |
| Engagement:Protest:Politics | 2.3180434 |
df_effects_activist <- get_all_predictions(
fit_activist$fit_updated[[19]], 'Activist support',
variables = c('Engagement', 'Protest', 'Politics')
)
p8 <- create_plot(
df, 'ActivistSupport', 'Activist support',
political = TRUE, df_effects = df_effects_activist
)
p8
ggsave('figures/FigureS5.pdf', p8, width = 9, height = 7)
inclusion_bfs_activist_support <- get_inclusion_bayes_factors(fit_activist)
kable(
log(inclusion_bfs_activist_support),
col.names = c('Effect', 'Log Inclusion Bayes factor')
)
| Effect | Log Inclusion Bayes factor |
|---|---|
| Engagement | -0.3266687 |
| Protest | 55.1913273 |
| Politics | 376.0767910 |
| Engagement:Protest | 0.7485845 |
| Engagement:Politics | -3.6072906 |
| Protest:Politics | -2.3995359 |
| Engagement:Protest:Politics | -0.3930084 |
df_effects_radical <- get_all_predictions(
fit_radical$fit_updated[[19]], 'Perceived radicalness',
variables = c('Engagement', 'Protest', 'Politics')
)
p9 <- create_plot(
df, 'Radical', 'Perceived radicalness', 'of the protest',
political = TRUE, df_effects = df_effects_radical,
lo = 'Not at all radical', hi = 'Extremely radical'
)
p9
ggsave('figures/FigureS6.pdf', p9, width = 9, height = 7)
inclusion_bfs_radical <- get_inclusion_bayes_factors(fit_radical)
kable(
log(inclusion_bfs_radical),
col.names = c('Effect', 'Log Inclusion Bayes factor')
)
| Effect | Log Inclusion Bayes factor |
|---|---|
| Engagement | 1.701844 |
| Protest | 205.525590 |
| Politics | 108.389725 |
| Engagement:Protest | -1.064140 |
| Engagement:Politics | -3.230613 |
| Protest:Politics | -2.129755 |
| Engagement:Protest:Politics | -1.442445 |
# Is small enough to be included, but also quick enough to be estimated once
if (!file.exists('models/fit_source_politics.RDS')) {
fit_source_full <- brm(
formula = 'Y ~ Engagement + Protest + Politics',
data = df_source, family = brms::cumulative('probit'),
prior = c(
prior(student_t(1, 0, 0.15), class = b),
prior(student_t(3, 0, 2.5), class = Intercept)
), cores = 4, chains = 4, iter = 2000, save_pars = save_pars(all = TRUE)
)
saveRDS(fit_source_full, 'models/fit_source_politics.RDS')
} else {
fit_source_full <- readRDS('models/fit_source_politics.RDS')
}
df_effects_source <- get_all_predictions(
fit_source_full, 'Source credibility',
variables = c('Engagement', 'Protest', 'Politics')
)
p10 <- create_plot(
df %>% filter(Engagement != 'Control'),
'SourceCredibility',
'Source credibility',
'Trust Dr. Fraser',
political = TRUE, lo = 'Not at all', hi = 'A great deal',
df_effects = df_effects_source
)
p10
ggsave('figures/FigureS7.pdf', p10, width = 9, height = 7)
inclusion_bfs_source <- get_inclusion_bayes_factors(fit_source)
kable(
log(inclusion_bfs_source),
col.names = c('Effect', 'Log Inclusion Bayes factor')
)
| Effect | Log Inclusion Bayes factor |
|---|---|
| Engagement | 0.2557725 |
| Protest | -0.8414234 |
| Politics | 186.9059814 |
| Engagement:Protest | -1.4305147 |
| Engagement:Politics | -2.1577978 |
| Protest:Politics | -1.0520638 |
| Engagement:Protest:Politics | -0.6796638 |
df_effects_science <- get_all_predictions(
fit_science$fit_updated[[19]], 'Science credibility',
variables = c('Engagement', 'Protest', 'Politics')
)
p11 <- create_plot(
df, 'ScienceCredibility',
'Science credibility',
'Trust environmental scientists',
political = TRUE, df_effects = df_effects_science,
lo = 'Not at all', hi = 'A great deal'
)
p11
ggsave('figures/FigureS8.pdf', p11, width = 9, height = 7)
inclusion_bfs_science <- get_inclusion_bayes_factors(fit_science)
kable(
log(inclusion_bfs_science),
col.names = c('Effect', 'Log Inclusion Bayes factor')
)
| Effect | Log Inclusion Bayes factor |
|---|---|
| Engagement | -3.346037 |
| Protest | -1.536171 |
| Politics | 421.843054 |
| Engagement:Protest | -2.281769 |
| Engagement:Politics | -3.599588 |
| Protest:Politics | -1.760329 |
| Engagement:Protest:Politics | -1.739139 |
fit_donation <- anovaBF(
Donation ~ Protest + Engagement + Politics, data = df,
rscaleEffects = c('Engagement' = 0.15, 'Protest' = 0.15, 'Politics' = 0.15)
)
if (!file.exists('models/df_effects_donation.csv')) {
df_preds_donation <- get_df_preds_donation(fit_donation[18], df)
df_effects_donation <- df_preds_donation %>%
group_by(Protest, Engagement, Politics) %>%
summarize(
estimate = mean(ypred),
conf.low = quantile(ypred, 0.025),
conf.high = quantile(ypred, 0.975),
)
write.csv(df_effects_donation, 'models/df_effects_donation.csv', row.names = FALSE)
} else {
df_effects_donation <- read.csv('models/df_effects_donation.csv')
}
df_effects_donation <- df_effects_donation %>%
mutate(
Engagement = ifelse(Engagement == 'Support', 'Endorsed', Engagement),
Engagement = factor(Engagement, labels = c('Control', 'Endorsed', 'Engaged'))
)
cols <- c('#08519c', 'gray46', '#a50f15')
pdon <- ggplot(
df_effects_donation, aes(x = Engagement, y = estimate, color = Politics)) +
geom_errorbar(
aes(ymin = conf.low, ymax = conf.high),
position = position_dodge(width = 0.70), width = 0.50, size = 1
) +
geom_point(size = 3, position = position_dodge(width = 0.70)) +
theme_minimal() +
labs(
x = '',
y = 'Predicted values',
title = 'Donations across conditions'
) +
scale_color_manual(values = cols) +
scale_y_continuous(limits = c(10, 40)) +
facet_wrap(~ Protest) +
theme(
legend.position = 'top',
legend.title = element_blank(),
strip.text = element_text(size = 12),
panel.grid = element_blank(),
plot.title = element_text(size = 14, hjust = 0.50)
) +
guides(fill = guide_legend(reverse = TRUE))
pdon
inclusion_bfs_donation <- c(
# Main effects
get_bf_inclusion_donation(fit_donation, 'Engagement')[1],
get_bf_inclusion_donation(fit_donation, 'Protest')[1],
get_bf_inclusion_donation(fit_donation, 'Politics')[1],
# Two-way interaction effects
get_bf_inclusion_donation(fit_donation, 'Engagement', varname2 = 'Protest')[2],
get_bf_inclusion_donation(fit_donation, 'Engagement', varname2 = 'Politics')[2],
get_bf_inclusion_donation(fit_donation, 'Protest', varname2 = 'Politics')[2],
# Tree-way interaction effects
get_bf_inclusion_donation(fit_donation, 'Engagement')[3]
)
names(inclusion_bfs_donation) <- names(inclusion_bfs_science)
kable(
log(inclusion_bfs_donation),
col.names = c('Effect', 'Log Inclusion Bayes factor')
)
| Effect | Log Inclusion Bayes factor |
|---|---|
| Engagement | -2.390265 |
| Protest | 1.545690 |
| Politics | 29.764425 |
| Engagement:Protest | -3.884956 |
| Engagement:Politics | -5.542848 |
| Protest:Politics | -4.276794 |
| Engagement:Protest:Politics | -3.976229 |
set.seed(1)
df_preds <- bind_rows(
get_all_predictions(fit_policy$fit_updated[[19]], outcome = 'Policy support'),
get_all_predictions(fit_activist$fit_updated[[19]], outcome = 'Activist support'),
get_all_predictions(fit_radical$fit_updated[[19]], outcome = 'Perceived radicalness'),
get_all_predictions(fit_science$fit_updated[[19]], outcome = 'Science credibility'),
get_all_predictions(fit_source_full, outcome = 'Source credibility'),
df_effects_donation %>% mutate(outcome = 'Donation')
) %>%
mutate(
Politics = factor(Politics, levels = c('Democrat', 'Independent', 'Republican')),
outcome = factor(
outcome,
levels = c(
'Policy support', 'Activist support',
'Perceived radicalness', 'Source credibility',
'Science credibility', 'Donation'
)
),
panel = c(rep(1:4, each = 18), rep(5, 12), rep(16, 18))
)
cols <- c('#08519c', 'gray46', '#a50f15')
ppost <- ggplot(
df_preds, aes(x = Engagement, y = estimate, col = Politics, shape = Protest)
) +
geom_errorbar(
aes(ymin = conf.low, ymax = conf.high),
position = position_dodge(width = 0.70), width = 0.50, size = 1
) +
geom_point(size = 3, position = position_dodge(width = 0.70)) +
theme_minimal() +
labs(
x = '',
y = 'Latent mean',
title = 'Latent means across outcomes and conditions'
) +
scale_color_manual(values = cols) +
facet_wrap(~ outcome, ncol = 2, scales = 'free_y') +
theme(
plot.title = element_text(size = 14, hjust = 0.50),
strip.text = element_text(size = 10),
legend.title = element_blank(),
legend.position = 'top',
legend.box = 'vertical',
legend.spacing.y = unit(-0.25, 'cm')
)
ggsave('figures/Figure2.pdf', ppost, width = 9, height = 9)
ppost
varnames <- c(
'Policy support',
'Activist support',
'Perceived radicalness',
'Source credibility',
'Science credibility',
'Donation'
)
inclusion_bfs <- cbind(
inclusion_bfs_policy_support,
inclusion_bfs_activist_support,
inclusion_bfs_radical,
inclusion_bfs_source,
inclusion_bfs_science,
inclusion_bfs_donation
)
format_scientific <- function(x) {
if (is.numeric(x) && x > 1000) {
formatted <- formatC(x, format = "e", digits = 2)
# Replace 'e+' with ' x 10^' for the desired format
formatted <- gsub("e\\+0*", " x 10^", formatted)
formatted <- gsub("e\\-", " x 10^-", formatted)
return(formatted)
}
return(as.character(x))
}
colnames(inclusion_bfs) <- varnames
df_inclusion <- data.frame(round(inclusion_bfs, 3)) %>% tibble::rownames_to_column()
colnames(df_inclusion) <- c('Factor', str_replace_all(colnames(df_inclusion)[-1], '\\.', ' '))
# Apply the formatting function to the data frame
df_inclusion[, -1] <- as.data.frame(lapply(df_inclusion[, -1], function(column) {
sapply(column, format_scientific)
}))
kable(df_inclusion)
| Factor | Policy support | Activist support | Perceived radicalness | Source credibility | Science credibility | Donation |
|---|---|---|---|---|---|---|
| Engagement | 0.036 | 0.721 | 5.484 | 1.291 | 0.035 | 0.092 |
| Protest | 0.369 | 9.32 x 10^23 | 1.81 x 10^89 | 0.431 | 0.215 | 4.691 |
| Politics | 3.65 x 10^158 | 2.13 x 10^163 | 1.18 x 10^47 | 1.49 x 10^81 | 1.60 x 10^183 | 8.44 x 10^12 |
| Engagement:Protest | 0.312 | 2.114 | 0.345 | 0.239 | 0.102 | 0.021 |
| Engagement:Politics | 0.056 | 0.027 | 0.04 | 0.116 | 0.027 | 0.004 |
| Protest:Politics | 0.108 | 0.091 | 0.119 | 0.349 | 0.172 | 0.014 |
| Engagement:Protest:Politics | 10.156 | 0.675 | 0.236 | 0.507 | 0.176 | 0.019 |
prior_widths <- seq(0.10, 0.50, 0.02) / 2
if (!file.exists('results/sensitivity_confirmatory.RDS')) {
# Policy support
h1_sens <- run_sensitivity(
df %>% mutate(Y = PolicySupport),
list('null' = fit_all[[1]], 'full' = fit_all[[3]]),
prior_widths
)
h4_sens <- run_sensitivity(
df %>% mutate(Y = PolicySupport),
list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
prior_widths, includes_ordered = TRUE,
groups = 3, pars = 'Engagement', type = 'increasing'
)
# Activist support
h2_sens <- run_sensitivity(
df %>% mutate(Y = ActivistSupport),
list('null' = fit_all[[1]], 'full' = fit_all[[3]]),
prior_widths
)
h5_sens <- run_sensitivity(
df %>% mutate(Y = ActivistSupport),
list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
prior_widths,
includes_ordered = TRUE,
groups = 3, pars = 'Engagement', type = 'increasing'
)
# Perceived radicalness
h3_sens <- run_sensitivity(
df %>% mutate(Y = Radical),
list('null' = fit_all[[1]], 'full' = fit_all[[3]]),
prior_widths
)
h6_sens <- run_sensitivity(
df_march,
list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
prior_widths,
groups = 3, pars = 'Engagement', type = 'decreasing'
)
h7_sens <- run_sensitivity(
df_civil,
list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
prior_widths,
groups = 3, pars = 'Engagement', type = 'decreasing'
)
# Source credibility
h8_sens <- run_sensitivity(
df_source,
list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
prior_widths
)
# Science credibility
h9_sens <- run_sensitivity(
df %>% mutate(Y = ScienceCredibility),
list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
prior_widths
)
sens_all <- list(
'h1_sens' = h1_sens,
'h2_sens' = h2_sens,
'h3_sens' = h3_sens,
'h4_sens' = h4_sens,
'h5_sens' = h5_sens,
'h6_sens' = h6_sens,
'h7_sens' = h7_sens,
'h8_sens' = h8_sens,
'h9_sens' = h9_sens
)
saveRDS(sens_all, 'results/sensitivity_confirmatory.RDS')
} else {
sens_all <- readRDS('results/sensitivity_confirmatory.RDS')
}
get_bf01 <- function(bfmat) bfmat[1, 2]
sapply(sens_all$h1_sens, get_bf01)
## [1] 3.086449 3.548767 4.045935 4.481988 5.049997 5.524717 6.076725
## [8] 6.485455 7.023044 7.590662 8.184941 8.699794 9.235297 9.596410
## [15] 10.182637 10.749493 11.251794 11.782758 12.399700 12.952234 13.457980
bfs <- do.call('c', lapply(names(sens_all), function(name) {
sapply(sens_all[[name]], get_bf01)
}))
df_sens <- data.frame(
type = rep(paste0('H[', seq(9), ']'), each = length(prior_widths)),
prior_width = rep(prior_widths, 9),
bfs = bfs
)
custom_label <- function(x) {
ifelse(
abs(x) > 500 | abs(x) < 0.002,
formatC(x, format = 'e', digits = 1),
format(x, scientific = FALSE)
)
}
psens <- ggplot(df_sens, aes(x = prior_width * 2, y = bfs)) +
geom_line(linewidth = 1) +
facet_wrap(~ type, scales = 'free_y', labeller = label_parsed) +
labs(
y = expression(BF[0*1]),
x = 'Prior width',
title = 'Bayes factor sensitivity analysis'
) +
scale_y_continuous(labels = custom_label) +
theme_minimal() +
theme(
legend.position = 'top',
legend.title = element_blank(),
strip.text = element_text(size = 12),
plot.title = element_text(size = 14, hjust = 0.50)
)
psens
ggsave('figures/FigureS2.pdf', psens, width = 9, height = 9)
prior_widths <- seq(0.10, 0.50, 0.04) / 2
policy_sens <- run_sensitivity_inclusion(
df %>% mutate(Y = PolicySupport),
fit_all,
prior_widths,
filename = 'results/policy_sens_incl.RDS'
)
activist_sens <- run_sensitivity_inclusion(
df %>% mutate(Y = ActivistSupport),
fit_all,
prior_widths,
filename = 'results/activist_sens_incl.RDS'
)
radical_sens <- run_sensitivity_inclusion(
df %>% mutate(Y = Radical),
fit_all,
prior_widths,
filename = 'results/radical_sens_incl.RDS'
)
science_sens <- run_sensitivity_inclusion(
df %>% mutate(Y = ScienceCredibility),
fit_all,
prior_widths,
filename = 'results/science_sens_incl.RDS'
)
source_sens <- run_sensitivity_inclusion(
df %>% mutate(Y = SourceCredibility),
fit_all,
prior_widths,
filename = 'results/source_sens_incl.RDS'
)
# For donation behaviour
donation_sens <- run_sensitivity_inclusion_donation(
df,
prior_widths,
filename = 'results/source_donation_incl.RDS'
)
all_sens <- rbind(
data.frame(policy_sens) %>% mutate(outcome = 'Policy support'),
data.frame(activist_sens) %>% mutate(outcome = 'Activist support'),
data.frame(radical_sens) %>% mutate(outcome = 'Perceived radicalness'),
data.frame(science_sens) %>% mutate(outcome = 'Science credibility'),
data.frame(source_sens) %>% mutate(outcome = 'Source credibility'),
data.frame(donation_sens) %>% mutate(outcome = 'Donation')
)
df_sens_all <- all_sens %>%
pivot_longer(cols = -c(prior_widths, outcome)) %>%
mutate(
name = gsub('\\.', ' x ', name),
name = factor(
name,
levels = c(
'Engagement', 'Protest', 'Politics',
'Engagement x Protest', 'Engagement x Politics',
'Protest x Politics', 'Engagement x Protest x Politics'
)
),
outcome = factor(
outcome,
levels = c(
'Policy support', 'Activist support',
'Perceived radicalness', 'Source credibility',
'Science credibility', 'Donation'
)
)
)
cols <- c('#D55E00', '#A65628', '#029E73', '#56B4E9', '#0072B2', '#CC79A7')
psens2 <- ggplot(df_sens_all, aes(x = prior_widths * 2, y = log(value), color = outcome)) +
geom_line(linewidth = 1) +
facet_wrap(~ name, scales = 'free_y') +
labs(
y = expression('Log Inclusion' ~ BF[1*0]),
x = 'Prior width',
title = 'Sensitivity analysis for exploratory inclusion Bayes factors'
) +
# scale_y_continuous(labels = custom_label) +
scale_color_manual(values = cols) +
theme_minimal() +
theme(
legend.position = 'top',
legend.title = element_blank(),
strip.text = element_text(size = 12),
plot.title = element_text(size = 14, hjust = 0.50),
legend.direction = "horizontal",
legend.box = "horizontal",
legend.box.just = "center"
) +
guides(colour = guide_legend(nrow = 2))
psens2
ggsave('figures/FigureS3.pdf', psens2, width = 9, height = 9)
df_cors <- df %>%
select(
PolicySupport, ActivistSupport, SourceCredibility, ScienceCredibility,
Radical, Donation, Gender, Age.x, Income, EducationLevel, PoliticalAffiliation
) %>%
mutate_all(as.numeric)
cormat <- cor(df_cors, method = 'kendall', use = 'pairwise.complete.obs')
diag(cormat) <- NA
rownames(cormat) <- colnames(cormat) <- c(
'Policy support', 'Activist support', 'Source credibility', 'Science credibility',
'Perceived radicalness', 'Donation', 'Gender', 'Age', 'Income', 'Educational level',
'Political affiliation'
)
pdf('figures/FigureS9.pdf', width = 8, height = 8)
corrplot(
cormat, method = 'color', number.cex = 0.80,
addCoef.col = 'black', type = 'upper',
tl.cex = 0.80, addrect = 20, tl.col = 'black',
na.label = ' ',
is.corr = TRUE,
number.digits = 2, col.lim = c(-0.80, 0.80)
)
dev.off()
## quartz_off_screen
## 2
corrplot(
cormat, method = 'color', number.cex = 0.80,
addCoef.col = 'black', type = 'upper',
tl.cex = 0.80, addrect = 20, tl.col = 'black',
na.label = ' ',
is.corr = TRUE,
number.digits = 2, col.lim = c(-0.80, 0.80)
)
sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.3
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Amsterdam
## tzcode source: internal
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2
## [4] stringr_1.5.1 marginaleffects_0.25.0 BayesFactor_0.9.12-4.7
## [7] Matrix_1.6-5 coda_0.19-4.1 gridExtra_2.3
## [10] corrplot_0.92 ggplot2_3.5.1 scales_1.3.0
## [13] readr_2.1.5 tidyr_1.3.1 knitr_1.45
## [16] dplyr_1.1.4 brms_2.21.0 Rcpp_1.0.12
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.1 farver_2.1.1 Rmpfr_0.9-5
## [4] loo_2.7.0 fastmap_1.1.1 tensorA_0.36.2.1
## [7] digest_0.6.35 estimability_1.5.1 lifecycle_1.0.4
## [10] StanHeaders_2.32.6 magrittr_2.0.3 posterior_1.5.0
## [13] compiler_4.3.3 rlang_1.1.3 sass_0.4.9
## [16] tools_4.3.3 utf8_1.2.4 yaml_2.3.8
## [19] collapse_2.1.0 data.table_1.17.0 labeling_0.4.3
## [22] bridgesampling_1.1-2 bit_4.0.5 pkgbuild_1.4.4
## [25] curl_5.2.1 abind_1.4-5 withr_3.0.0
## [28] purrr_1.0.2 grid_4.3.3 stats4_4.3.3
## [31] fansi_1.0.6 xtable_1.8-4 colorspace_2.1-0
## [34] inline_0.3.19 emmeans_1.10.2 insight_1.1.0
## [37] cli_3.6.2 mvtnorm_1.2-4 rmarkdown_2.27
## [40] crayon_1.5.2 ragg_1.3.0 generics_0.1.3
## [43] RcppParallel_5.1.7 rstudioapi_0.16.0 tzdb_0.4.0
## [46] pbapply_1.7-2 cachem_1.0.8 rstan_2.32.6
## [49] bayesplot_1.11.1 matrixStats_1.3.0 vctrs_0.6.5
## [52] V8_4.4.2 jsonlite_1.8.8 hms_1.1.3
## [55] bit64_4.0.5 systemfonts_1.0.6 jquerylib_0.1.4
## [58] glue_1.7.0 codetools_0.2-19 distributional_0.4.0
## [61] stringi_1.8.3 gtable_0.3.5 QuickJSR_1.1.3
## [64] gmp_0.7-4 munsell_0.5.1 tibble_3.2.1
## [67] pillar_1.9.0 htmltools_0.5.8.1 Brobdingnag_1.2-9
## [70] R6_2.5.1 textshaping_0.3.7 vroom_1.6.5
## [73] evaluate_0.23 lattice_0.22-5 highr_0.10
## [76] backports_1.4.1 bslib_0.7.0 rstantools_2.4.0
## [79] MatrixModels_0.5-3 nlme_3.1-164 checkmate_2.3.1
## [82] xfun_0.43 pkgconfig_2.0.3